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^ i 1. Overview 



Algebraic manipulation on computer is a tool that has been developed quite soon, about 



m 

^ . one decade after the birth of computers, the first examples dating back to the end of 
the fifties of the last century. General purpose packages began to be developed during 
the sixties, and include, e.g.. Reduce (1968), Macsyma (1978), muMath (1980), Maple 
; (1984), Scratchpad (1984), Derive (1988), Mathematica (1988), Pari/GP (1990) and 

^ ■ Singular (1997) (the dates refer to the first release). However, most of the facilities of 
■ ■ ■ these general purpose manipulators are simply ignored when dealing with perturbation 

methods in Celestial Mechanics. For this reason, the job of developing specially devised 
manipulation tools has been undertaken by many people, resulting in packages that have 
limited capabilities, but are definitely more effective in practical applications. Producing 
a list of these packages is a hard task, mainly because most of them are not publicly 
available. A list of "old time" packages may be found in Henrard [18] and Laskar [21]. 
In recent times a manipulator developed by J. Laskar and M. Gastineau has become 
quite known. 

Finding references to the methods implemented in specially devised packages is as 
difficult as giving a list. We know only a few papers by Broucke and Garthwaite [3], 
Broucke [4], Rom [24], Henrard [17] and [18], Laskar [21], Jorba [19] and Biscani [2]. 
A complete account of the existing literature on the subject goes beyond the limits of 
the present note. The present work introduces some ideas that have been used by the 
authors in order to implement a package named Xgouot,. 
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As a matter of fact, most of the algebraic manipulation used in Celestial Mechanics 
makes use of the so called "Poisson series", namely series with a general term of the 
form 

• . . . • x^" . {.ki^pi + . . . + kmVm) , 
sm 

(with obvious meaning of the symbols). Thus, a very minimal set of operations is re- 
quired, namely sums, products and derivatives of polynomials and/or trigonometric 
polynomials. Traditionally, also the operation of inversion of functions, usually made 
again via series expansion, was required. However, the expansion methods based on Lie 
series and Lie transforms typically get rid of the latter operation (see, e.g., [16]). 

Writing a program doing algebraic manipulation on series of the type above leads 
one to be confronted with a main question, namely how to represent a polynomial, 
trigonometric polynomial or Poisson series on a computer. The papers quoted above 
actually deal with this problem, suggesting some methods. In these lectures we provide 
an approach to this problem, followed by a few examples of applications. 

In sect. 2 we include a brief discussion about the construction of normal form for a 
Hamiltonian system in the neighborhood of an elliptic equilibrium. We do not attempt 
to give a complete discussion, since it is available in many papers. We rather try to 
orient the reader's attention on the problem of representing perturbation series. 

In sect. 3-7 we introduce a method which turns out to be quite useful for the 
representation of a function as an array of coefficients. The basic idea has been suggested 
to one of the authors by the paper of Gustavson [14] (who, however, just mentions that 
he used an indexing method, without giving any detail about its implementation). One 
introduces an indexing function which transforms an array of exponents in a polynomial 
(or trigonometric polynomial) in a single index within an array. The general scheme 
is described in sect. 3. The basics behind the construction of an indexing function 
are described in sect. 4. The details concerning the representation of polynomials and 
trigonometric polynomials are reported in sects. 5 and 6, respectively. In sect. 7 we 
include some hints about the case of sparse series, that may be handled by combining 
the indexing functions above with a tree representation. Finally, sect. 8 is devoted to 
three applications, by giving a short account of the contents of published papers. 



2. A common problem in perturbation theory 

A typical application of computer algebra is concerned with the construction of first 
integrals or of a normal form for a Hamiltonian system. A nontrivial example, which 
however may be considered as a good starting point, is the calculation of a normal 
form for the celebrated model of Henon and Heiles [15], which has been done by Gus- 
tavson [14]. Some results on this model are reported in sect. 8. 

We assume that the reader is not completely unfamiliar with the concept of normal 
form for a (possibly Hamiltonian) system of differential equations. Thus, let us briefly 
illustrate the problem by concentrating our attention on the algorithmic aspect and by 
explaining how algebraic manipulation may be introduced. 
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2.1 Computation of a normal form 

Let us consider a canonical system of differential equations in the neighborhood of an 
elliptic equilibrium. The Hamiltonian may typically be given the form of a power series 
expansion 



n 

UJ-i 



(1) H(x,y) = Ho(x,y) + Hi(x,y) + ... , ffola;, y) = J] ^(xj + , 

3 = 1 

where Hs{x, y) for s > 1 is a homogeneous polynomial of degree s + 2 in the canonical 
variables (x, y) G M^". Here a; e is the vector of the frequencies, that are assumed 
to be all different from zero. 

In such a case the system is said to be in Birkhoff normal form in case the Hamil- 
tonian takes the form 

(2) H{x,y) = Ho{x,y) + Zi{x,y) + Z2{x,y) + . . . with L//„Zs = , 

where Lhq- = {Hq, ■} is the Lie derivative with respect to the flow of Hq, actually the 
Poisson bracket with Hq. 

The concept of Birkhoff normal form is better understood if one assumes also that 
the frequencies are non resonant, i.e., if 

{k, u)^0 for all e , A; 7^ , 

where (/c, u) = kjuij. For, in this case the functions Zs{x, y) turn out to be actually 
function only of the n actions of the system, namely of the quantities 

x] + y] . 
Ij= 2 ' 3 = l,...,n. 

It is immediate to remark that /i , . . . , are independent first integrals for the Hamil- 
tonian, an that they are also in involution, so that, by Liouville's theorem, the system 
turns out to be integrable. The definition of normal form given in (2) is more general, 
since it includes also the case of resonant frequencies. 

The calculation of the normal form may be performed using the Lie transform 
method, which turns out to be quite effective. We give here the algorithm without proof. 
A complete description may be found, e.g., in [9], and the description of a program 
implementing the method via computer algebra is given in [10]. The corresponding 
FORTRAN program is available from the CPC library. 

The Lie transform is defined as follows. Let a generating sequence xi(a;, y), 
X2{x, y), ... he given, and define the operator 



(3) r^ = E^^ 

where the sequence Eq, Ei, . . . of operators is recursively defined as 

s 

(4) £^0 = 1, Es = J2^-L^^Es 



S ' 
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This is a linear operator that is invertible and satisfies the interesting properties 

(5) T^{f,g) = {TJ,T^g) , T^{f ■ g) = ■ T^g . 

Let now Z{x, y) = Hq{x, y) + Zi{x, y) + Z2{x,y) + . . . be a function such that 

(6) T^Z = H , 

where H is our original Hamiltonian, and let Z possess a first integral $, i.e., a function 
satisfying {Z, $} = 0. Then one has also 

T^{Z, $} ^ {T^Z, T^$} = {H, T^$} = , 

which means that if $ is a Grst integral for Z then is a iirst integral for H. 

The question now is: can we find a generating sequence xi, Xii - ■ ■ such that the 
function Z satisfying (6 ) is in Birkhoff normal form ? 

The answer to this question is in the positive, and the generating sequence may be 
calculated via an explicit algorithm that can be effectively implemented via computer 
algebra. We include here the algorithm, referring to, e.g., [9] for a complete deduction. 
Here we want only to stress that all operations that are required may be actually 
implemented on a computer. 

The generating sequence is determined by solving for x ctnd Z the equations 

(7) Zs - LhoXs = Hs + Qs , s>l, 

where Qs is a known homogeneous polynomial of degree s + 2 given by = and 

s-l 

Qs = -J2iEjZs-j + -{Xj, Es-jHo}) , s>l. 
i=i ^ 

In order to solve (7) it is convenient to introduce complex variables ^, rj via the canonical 
transformation 

1 i 

which transforms Hq — i ujjCjVj- these variables the operator Lh^ takes a diagonal 
form, since 

where we have used the multi-index notation = • . . . • and similarly for rj. 
Thus, writing the r.h.s. of (7) as a sum of monomials Cj^k^^r]^ the most direct form of 
the solution is found by including in Z all monomials with {k — j, oj) = 0, and adding 
i{k-'j uj) ^'' ^° ^'-"^ monomials with {k — j,u)) ^ 0. This is the usual way of 
constructing a normal form for the system (1). 

Let us now examine in some more detail the algebraic aspect. With some patience 
one can verify that (7) involves only homogeneous polynomials of degree s + 2. Thus, one 
should be able to manipulate this kind of functions. Moreover, a careful examination 
of the algorithm shows that there are just elementary algebraic operations that are 
required, namely: 

(i) sums and multiplication by scalar quantities; 
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(ii) Poisson brackets, which actuaUy require derivatives of monomials, sums and prod- 
ucts; 

(iii) hnear substitution of variables, which may still be reduced to calculation of sums 
and products without affecting the degree of the polynomial; 

(iv) solving equation (7), which just requires a division of coefficients. 

These remarks should convince the reader that implementing the calculation of the 
normal form via algebraic manipulation on a computer is just matter of being able of 
representing homogeneous polynomials in many variables and performing on them a few 
elementary operations, such as sum, product and derivative. 

2.2 A few elementary considerations 

In order to have an even better understanding the reader may want to consider the 
elementary problem of representing polynomials in one single variable. We usually write 
such a polynomial of degree s (non homogeneous, in this case) as 

f{x) = ao + aix + . . . + agX^ . 

A machine representation is easily implemented by storing the coefficients ao, oi, . . ., On 
as a one-dimensional array of floating point quantities, either real or complex. E.g., in 
FORTRAN language one can represent a polynomial of degree 100 by just saying, e.g., 
DIMENSION F(lOl) and storing the coefficient aj as F(j+1) (here we do not use the 
extension of FORTRAN that allows using zero or even negative indices for an array). 
Similarly in a language like C one just says, e.g., double f [101] and stores aj as f [j] . 

The operation of sum is a very elementary one: if /, g are two polynomials and 
the coefficients are stored in the arrays f , g (in C language) then the sum h is the 
array h with elements hCj] = f [j] + g[j]. The derivative of / is the array fp with 
elements fp[j] = (j+l)*f [j+1] . In a similar way one can calculate the product, by 
just translating in a programming language the operations that are usually performed 
by hand. 

The case of polynomials in two variables is just a bit more difficult. A homogeneous 
polynomial of degree s is usually written as 

f{x, y) = as,ox^ + as-i,ix^y + . . . + ao,s2/* • 

The naive (not recommended) representation would use an array with two indices (a 
matrix), by saying, e.g., DIMENSION F(101,101) and storing the coefficient aj^k as 
F(j+l,k+l). Then the algebra is just a straightforward modification with respect to 
the one-dimensional case. 

Such a representation is not recommended for at least two reasons. The first one 
is that arrays with arbitrary dimension are difficult to use, or even not allowed, in 
programming languages. The second and more conclusive reason is that such a method 
turns out to be very effective in wasting memory space. E.g., in the two dimensional 
case a polynomial of degree up to .s requires a matrix with {s + 1)^ elements, while only 
(s -|- l)(s + 2)/2 are actually used. Things go much worse in higher dimension, as one 
easily realizes. 



6 



The arguments above should have convinced the reader that an effective method 
of representing polynomials is a basic tool in order to perform computer algebra for 
problems like the calculation of normal form. Once such a method is available, the rest 
is essentially known algebra, that needs to be translated in a computer language. 

The problem for Poisson series is a similar one, as the reader can easily imagine. 
The following sections contains a detailed discussion of indexing methods particularly 
devised for polynomials and for Poisson series. The underlying idea is to represent the 
coeflBcients as a one-dimensional array by suitably packing them in an effective manner, 
so as to avoid wasting of space. 

3. General scheme 

The aim of this section is to illustrate how an appropriate algebraic structure may help 
in representing the particular classes of functions that usually appear in perturbation 
theory. We shall concentrate our attention only on polynomials and trigonometric poly- 
nomials, which are the simplest and most common cases. However, the reader will see 
that most of the arguments used here apply also to more general cases. 

3. 1 Polynomials and power series 

Let V denote the vector space of polynomials in the independent variables x = 
{xi, . . . , Xn) e M"". A basis for this vector space is the set {'Ufc(^)}fc6Z" , where 



In particular, we shall consider the subspaces Vs of V that contain all homogeneous 
polynomials of a given degree s > 0; the subspace Vo is the one-dimensional space of 
constants, and its basis is {1}. The relevant algebraic properties are the following: 

(i) every subspace Vs is closed with respect to sum and multiplication by a number, 
i.e., ilfeVs A geVs then f + g e Vs and af e Vs; 

(ii) the product of homogeneous polynomials is a homogeneous polynomial, i.e., if / G 
Vr A geVs then fg G Vr+s', 

(iii) the derivative with respect to one variable maps homogeneous polynomials into 
homogeneous polynomials, i.e., if / e then dx f € Vs-i', if s = then f — 0, 



These three properties are the basis for most of the algebraic manipulations that are 
commonly used in perturbation theory. 

A power series is represented as a sum of homogeneous polynomials. Of course, in 
practical calculations the series will be truncated at some order. Since every homoge- 
neous polynomial f gVs can be represented as 



(8) 



Uk{x) = x'' = x^^' - ...-x^ 



of course. 




\k\=s 



it is enough to store in a suitable manner the coefficients ft- A convenient way, par- 
ticularly effective when most of the coefficients are different from zero, is based on the 



7 



Table 1. Illustrating the function representation for power series. A memory block 
is assigned to the function f(x). The coefficient fk of Uk{x) is stored at the address 
resulting by adding the offset I{k) to the starting address of the memory block. 



/(0,0v,0) 



/(i,o,-,o) 



/(fcl, 



^2 ; • • • ;fcn ) 



^0 = I{k)^k = (0,0,..., 0) 
^ 1 = I(k) ^ A;= (1,0,...,0) 



I{k) k = {ki, k2,.. ., kn) 



usual lexicographic ordering of polynomials (to be pedantic, inverse lexicographic). E.g., 
a homogeneous polynomial of degree s in two variables is ordered as 

asfixl + as-i,ixl~^X2 + . . . + 00,5X2 . 

The idea is to use the position of a monomial in the lexicographic order as an 
index /(/ci, . . . , k^) in an array of coefficients. We call / and indexing function. Here we 
illustrate how to use it, deferring to sect. 5 the actual construction of the function. 

The method is illustrated in table 1. Let / be a power series, truncated at some 
finite order s. A memory block is assigned to /. The size of the block is easily determined 
as /((O, . . . , 0, s)). For, (0, . . . , 0, s) is the last vector of length s. The starting address 
of the block is assigned to the coefficient of 'U(o,o,...,o); the next address is assigned to 
the coefficient of 'U(i,o,...,o); because (1,0,..., 0) is the first vector of length 1, and so on. 
Therefore, the address assigned to the coefficient of tt(fci,...,fc„) is the starting address of 
the block incremented by /((fci, . . . , kn))- If / is a homogeneous polynomial of degree 
s the same scheme works fine with a few minor differences: the length of the block is 
/((O, . . . , 0, s)) — /((O, . . . , 0, s — 1)), the starting address of the block is associated to the 
coefficient of W(s,o,...,o)7 cind the coefficient of 'U(fci,...,fc„) is stored at the relative address 
/((/ci, . . . , kn)) — /((O, . . . , 0, s — 1)) . This avoids leaving an empty space at the top of 
the memory block. 

In view of the form above of the representation a function is identified with a set 
of pairs (/c, fk), where k G is the vector of the exponents, acting as the label of the 
elements of the basis, and fk is the numerical coefficient. Actually the vector k is not 
stored, since it is found via the index. The algebraic operations of sum, product and 
differentiation can be considered as operations on the latter set. 
(i) If /, 5r e Va then the operation of calculating the sum / + is represented as 



{k, fk 
{k, gk 

to be executed over all k such that \k 



^ {k,fk + 9k) , 
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(ii) If f EVr and g &Vs then the operation of calculating the product fg is represented 
as 

fk) ) ^ {k + k'Jkgk') , 
{k , gk')} 

to be executed over all k, k' such that \k\ = r and \k'\ = s. 

(iii) If f E Vs then the operation of differentiating / with respect to, e.g., xi is repre- 
sented as 

r0 for fci = , 

^ '^''"^{{k'^kM for /ci 7^0, 
where fc' = {ki — 1, ^2, • • • , ^n)- 

It is perhaps worthwhile to spend a few words about how to make the vector k to 
run over all its allowed values. In the case of sum, we do not really need it: since the 
indexes of both addends and of the result are the same, the operation can actually be 
performed no matter which k is involved: just check that the indexes are in the correct 
range. ^ In order to perform product and differentiation it is essential to know the values 
of k and k' . To this end, we can either use the inverse of the indexing function, or 
generate the whole sequence by using a function that gives the vector next to a given k. 

3.2 Fourier series 

Let us denote hy (p = {cpi . . . , (pn) G T"' the independent variables. The Fourier expan- 
sion of a real function on T'* takes the form 

(9) /(¥')= {akCOs{k,(p) + bksm{k,(p)) , 

where and bk are numerical coefficients. In this representation there is actually a lot 
of redundancy: in view of cos(— a) = coso: and sin(— a) = — sino; the modes —k and k 
can be arbitrarily interchanged. On the other hand, it seems that we actually need two 
different arrays for the sin and cos components, respectively. A straightforward way out 
is to use the exponential representation '^^.ake'^^''''^\ but a moment's thought leads us 
to the conclusion that the redundancy is not removed at all. However, we can at the 
same time remove the redundancy and reduce the representation to a single array by 
introducing a suitable basis {uk{ip)}keZ" • Let k e Z"; we shall say that k is even if 
the first non zero component of k is positive, and that k is odd if the first non zero 
component of k is negative. The null vector A; = is said to be even. Then we set 

rcos(A;, (fi) for k even , 
\sin(A;,<^) for A; odd . 

This makes the representation f{(fi) = Y2keZ" V'ktJ'ki'fi) unique and redundancy free. 
It may be convenient to remark that the notation for the sin function may create 

^ For a homogeneous polynomial of degree s the first vector is (s, 0, . . . , 0), and the last one 
is (0, ... , 0, s). The indexes of these two vectors are the limits of the indexes in the sum. 
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some confusion. Usually, working with one variable, we write sinip. The convention 
above means that we should rather write — sin(— (/p), which is correct, but a bit funny. 
This should be taken into account when, after having accurately programmed all the 
operations, we discover that our manipulator says, e.g., that ^ cos<^ = — sin(— (^). 

In view of the discussion in the previous section it should now be evident that a 
truncated Fourier expansion of a function f{ f) can easily be represented by storing the 
coefficient of Uk{(fi) at an appropriate memory address, as calculated by the indexing 
function I{k) of sect. 6. 

The considerations of the previous section can be easily extended to the problem of 
calculating the sum and/or product of two functions, and of differentiating a function. 
Let us identify any term of the Fourier expansion of the function / with the pair {k, fk)- 
Let us also introduce the functions odd(/c) and even(/c) as follows: if k is odd, then 
odd(/c) = k and even(/c) = —k; else odd(/c) = —k and even(/c) = k. That, is, force k to 
be odd or even, as needed, by possibly changing its sign. 

(i) Denoting by (k, fk) and {k,gk) the same Fourier components of two functions / 
and g, respectively, the sum is computed as 



(11) 



{k, fk) 
{k, Qk) 



^ {k, fk + 9k) 



(ii) Denoting by {k, fk) and {k',gk') any two terms in the Fourier expansion of the 
functions / and g, respectively, the product is computed as 



(12) 



(fc, fk) 
{k', gk>) 



eveii{k + /c'), y | even(A; — A;') '^'^^^ 



odd(A; + k'), Ih^^ u (^odd(/c - k') 



2 

fkgk 



odd{k + k'), ] u ( odd{k - k 



2 

/N fkOk' 



for k even , k' even , 
for k even , k' odd , 
for k odd , k' even , 



(^even{k + k'), j u (^eyen{k - k'), for k odd , k' odd . 

Remark that the product always produces two distinct terms, unless A; = or A;' = 0. 
(iii) Denoting by (k, fk) any term in the Fourier expansion of a function /, differentiation 
with respect to, e.g., is performed as 



(13) 



{kjk) ^ 



{—k, —kifk) for k even , 
{—k, kifk) for k odd . 



All these formulae follow from well known trigonometric identities. 
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4. Indexing functions 

The basic remark for constructing an index function is the foUowing. Suppose that we 
are given a countable set A. Suppose also that A is equipped with a relation of complete 
ordering, that we shall denote by the symbols -<, zii, y and >z. So, for any two elements 
a,b E A exactly one of the relations a -< h, a = h and 6 >- a is true. Suppose also that 
there is a minimal element in A., i.e., there is ao € ^ such that a )^ for all a e ^ 
such that a ^ ttQ. Then an index function I is naturally defined as 

(14) I{a) = 4^{heA:h~<a}. 

If ^ is a finite set containing N elements, then I{Ai) — {0, 1, . . . , AT — 1}. If ^ is an 
infinite (but countable) set, then I {A) = Z+, the set of non negative integers. For 
instance, the trivial case is ^ = Z_|_ equipped with the usual ordering relation. In such 
a case the indexing function is just the identity. 

Having defined the function /(a), we are interested in performing the following basic 
operations: 

(i) for a given a E A, find the index /(a); 

(ii) for a given a E A^ find the element next (or prior) to a, if it exists; 

(iii) for a given I G /(^), find /~^(/), i.e., the element a & A such that /(a) = /. 

The problem here is to implement an effective construction of the index for some 
particular subsets of Z"^ that we are interested in. In order to avoid confusions, we shall 
use the symbols -<, ^, >- and >z when dealing with an ordering relation in the subset 
of 17^ under consideration. The symbols <, <, > and > will always denote the usual 
ordering relation between integers. 

As a first elementary example, let us consider the case ^ = Z. The usual ordering 
relation < does not fulfill our requests, because there is no minimal element. However, 
we can construct a different ordering satisfying our requests as follows. 

Let k,k' G Z. We shall say that k' ^ k in case one of the following relations is true: 

(i) \k'\ < \k\ ; 

(ii) \k'\ = \k\ A k' >k. 

The resulting order is 0, 1, —1, 2, —2, . . . , so that is the minimal element. 
Constructing the indexing function in this case is easy. Indeed, we have 



2a — 1 for a > , 
—2a for a < . 
The inverse function is also easily constructed: 

■(/ + l)/2 fori odd, 
-1/2 for I even . 



(15) /(0) = 0, /(a) 



(16) I-\0)^0, /-^(0 = {| 



In the rest of this section we show how an indexing function can be constructed for 

two particularly interesting cases, namely polynomials and trigonometric polynomials. 
However, we stress that the procedure we are using is a quite general one, so it can be 
extended to other interesting situations. 
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5. The polynomial case 

Let us first take An = Z", i-e., integer vectors with non negative components; formally 

An = {k = {ki, . . . , fcn) e : fci > 0, . . . , fcn > 0} . 

The index n in An denotes the dimension of the space. This case is named "polynomial" 
because it occurs precisely in the representation of homogeneous polynomials, and so 
also in the Taylor expansion of a function of n variables: the integer vectors k e An 
represent all possible exponents. 

We shall denote hj \k\ = ki + . . . + kn the length (or norm) of the vector k G Z!J:. 
Furthermore, to a given vector k — (/ci,...,fcn) £ An we shall associate the vector 
t{k) e An-i (the tail of k) defined as t{k) = (/c2, • • • , kn)- This definition is meaningful 
only if n > 1, of course. 

5.1 Ordering relation 

Pick a fixed n, and consider the finite family of sets Ai = , . . . , An = 

Let k,k' G Am, with any 1 < m < n. We shall say that k' -< k in case one of the 
following conditions is true: 

(i) m > 1 A \k'\ < \k\ ; 

(a) m > 1 A \k'\ = \k\ A k[ > ki ; 
(Hi) m> 1 A \k'\ = \k\ A k[ = ki A t{k') -< t{k) . 

In table 2 the ordering resulting from this definition is illustrated for the cases m = 
2,3,4,5. 

If m = 1 then only (i) applies, and this ordering coincides with the natural one in 
Z-(.. For m> 1, if (i) and (ii) are both false, then (iii) means that one must decrease the 
dimension n by replacing k with its tail t{k), and retry the comparison. For this reason 
the ordering has been established for 1 < m < n. Eventually, one ends up with m = 1, 
to which only (i) applies. 

It is convenient to define Vn{k) as the set of the elements which precede k\ formally: 

Vn{k) = {k' eAn : k' -<k} . 

With the latter notation the indexing function is simply defined as I{k) = ^Vn{k). The 
following definitions are also useful. Pick a vector k G Am and define the sets Bn\k), 
Bn^\k) and Bn"\k) as the subsets of An satisfying (i), (ii) and (iii), respectively, in 
the ordering algorithm above. Formally: 

i3«(0) = Bl^'\0) = Bl^"\0) = B["\k) = B^^''\k) = ij) , 

Bi:Hk) = {k' e An : \k'\<\k\}, 
Bl^'^k) = {k' G An : = \k\ A k[ > ki} , 
B'^^'^k) = {k' G An : \k'\ = \k\ A k[ = ki A t{k') < t{k)} . 

The sets Bi^\k), B'n'Hk) and Bli''\k) are pairwise disjoint, and moreover 

B^\k) U B^^^\k) U Bl:^^\k) ^ Vnik) . 



Table 2. Ordering of integer vectors in for m = 2, 3, 4, 5. 
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(0, 


2, 


0, 


0, 
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0, 


1) 
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1, 


1, 


0, 
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1, 
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1, 
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1, 
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0, 
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0, 
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2, 
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0, 


2) 



This easily follows from the definition. 

5.2 Indexing function for polynomials 

Let k G An- In view of the definitions and of the properties above the index function, 
defined as in (14), turns out to be 

(18) /(O) = , I{k) = #4^) (k) + #4-) {k) + #4-^) {k) . 

Let us introduce the functions 

J{n,s) = e An ■■ \k\ = s} , 

N{n,s) = ^J{nJ) forn>l,s>0. 
These functions will be referred to in the following as J-table and N-table. 



13 



We claim that the indexing function can be recursively computed as 

m = , 

(20) f N(n, \k\ - 1) for ki = \k\ , 
^ ^ I{k) = J ' ' ' ' 

[N{n,\k\-1) + I{t{k)) foTki<\k\. 

The claim follows from 

(21) #B«(fc) = Ar(n,|A;|-l) ; 

for ki = \k\ , 

N{n-l,\k\ -ki-1) for ki < \k\ ; 



(22) #4"^/^) = 

(23) #B(:")(A;) = 



Iit{k)) for A;i = |A;| , 

litik)) - N{n - 1, \k\ -ki-1) for ki < \k\ . 

The equality (21) is a straightforward consequence of the definition of the A'^-table. 
The equality (22) follows from (17). Indeed, for |A;| = ki we have Bn^\k) = 0, and for 
\k\ > ki we have 

^ii'Kk) = U {k' eAn : k[=j A \t{k')\ = \k\-j} 
ki<j<\k\ 

[j {k' e An : k[ = \k\ - I A \t{k')\ = 1} . 

0<l<\k\-ki 

Coming to (23), first remark that 

B!^''\k) = {k' eAn:k[ = kiA \t{k')\ = \k\-ki A t{k') ■< t{k)} , 

so that 

#i3i"*)(/c) = #{A e An-i : |A| = |A;| - fci A A ^ t{k)} . 
Then, the equality follows by remarking that 

Vn-i{tik)) = {A e An-i ■■ |A| = \k\-ki A A ^ t{k)} U {A e An-i ■■ |A| < - ki} . 
Adding up all contributions (20) follows. 

5.3 Construction of the tables 

In view of (19) and (20) the indexing function is completely determined in explicit form 
by the J-table. We show now how to compute the J-table recursively. For n = 1 we 
have, trivially, J(l, s) = 1 for s > 0. For n> 1 use the elementary property 

{k e An ■■ \k\^ s} ^ [j {keAn ■■ ki = s-j A \tik) \ = j} . 

o<i<s 

Therefore 

J(l,s) = l , 

^^^"^ J{n,s) = ^J{n - for n > 1 . 

j=0 
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This also means that, according to (19), we have N{n, s) = J{n + 1, s) . 

By the way, one will recognize that the J-table is actually the table of binomial 
coefficients, being J{n,s) — ^"^^7^) • 

5.4 Inversion of the index function 

The problem is to find the vector k E corresponding to a given index I. 

For n = 1 we have = l, of course. Therefore, let us assume n > 1. We shall 

construct a recursive algorithm which calculates the inverse function by just showing 
how to determine ki and l(t{k)). 

(i) If / = 0, then /c = 0, and there is nothing else to do. 

(ii) If / > 0, find an integer s satisfying N{n, s — 1) < I < N{n, s). In view of (20) we 
have \k\ = s and I{t{k)) = I — N{n,s — 1). Hence, by the same method, we can 
determine |t(A;)|, and so also ki = s — \t{k)\. 

5.5 An example of implementation 

We include here an example of actual implementation of the indexing scheme for poly- 
nomials. This is part of a program for the calculation of first integrals that is fully 
described in [10]. The complete computer code is also available from the CPC program 
library. 

We should mention that the FORTRAN code included here has been written in 1976. 
Hence it may appear a little strange to programmers who are familiar with the nowadays 
compilers, since it does not use many features that are available in FORTRAN 90 or in 
the current versions of the compiler. It rather uses the standard of FORTRAN II, with 
the only exception of the statement PARAMETER that has been introduced later. 

The PARAMETERS included in the code allow the user to control the allocation of 
memory, and may be changed in order to adapt the program to different needs. 
NPMAX is the maximum number of degrees of freedom 
NORDMX is the maximal polynomial degree that will be used 

NBNl and NBN2 are calculated from the previous parameters, and are used in order to 
allocate the correct amount of memory for the table of binomial coefficients. Here are 
the statements: 

PARAMETER (NPMAX=3) 
PARAMETER (N0RDMX=40) 
PARAMETER (NBN2=2*NPMAX) 
PARAMETER (NBN1=N0RDMX+NBN2) 

As explained in the previous sections, the indexing function for polynomials uses 
the table of binomial coefficients. The table is stored in a common block named BINTAB 
so that it is available to all program modules. In the same block there are also some 
constants that are used by the indexing functions and are defined in the subroutine 
BINOM below. Here is the statement that must be included in every source module that 
uses these data: 

COMMON /BINTAB/ IBIN (NBNl , NBN2) , NPIUl , NMENl , NFAT , NBIN 
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Subroutine BINOM fills the data in the common block BINTAB. It must be called at 
the beginning of the execution, so that the constants become available. Forgetting this 
call will produce unpredictable results. The calling arguments are the following. 
NLIB : the number of polynomial variables. In the Hamiltonian case considered in the 
present notes it must be set as 2n, where n is the number of degrees of freedom. It must 
not exceed the value of the parameter NPMAX. 

NORD : the wanted order of calculation of the polynomials, which in our case is the 
maximal order of the normal form. It must not exceed the value of the parameter NFAT. 

The subroutine checks the limits on the calling arguments; if the limits are violated 
then the execution is terminated with an error message. The calculation of the part of 
the table of binomial coefficients that will be used is based on well known formulae. 
SUBROUTINE BINOM (NLIB, NORD) 

C 

C Compute the table of the binomial coefficients. 

C 

COMMON /BINTAB/ IBIN (NBNl , NBN2) , NPIUl , NMENl , NFAT , NBIN 

C 

NFAT=NORD+NLIB 
NBIN=NLIB 

IF (NFAT. GT. NBNl. OR. NBIN. GT.NBN2) GO TO 10 

NPIUl = NLIB+1 

NMENl = NLIB-1 

DO 1 1=1, NFAT 

IBIN(I,1) = I 

DO 1 K=2,NBIN 

IF(I-K) 2,3,4 

2 IBIN(I,K) = 
GO TO 1 

3 IBIN(I,K) = 1 
GO TO 1 

4 IBIN(I,K) = IBIN(I-1,K-1)+IBIN(I-1,K) 
1 CONTINUE 

RETURN 

10 WRITE (6, 1000) NFAT, NBIN 
STOP 

1000 FORMAT (//,5X,15HERR0R SUB BINOM, 2110,//) 

END 

Function INDICE implements the calculation of the indexing function for polynomi- 
als. The argument J is an integer array of dimension NLIB which contains the exponents 
of the monomial. It must contain non negative values with sum not exceeding the value 
NORD initially passed to the subroutine BINOM. These limits are not checked in order to 
avoid wasting time: note that this function may be called several millions of times in 
a program. The code actually implements the recursive formula (20) using iteration. 
Recall that recursion was not implemented in FORTRAN II. 
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FUNCTION INDICE(J,NLIB) 

C 

C Compute the relative address I corresponding to the 

C exponents J. 

C 

COMMON /BINTAB/ IBIN (NBNl , NBN2) , NPIUl , NMENl , NFAT , NBIN 
DIMENSION J(NLIB) 

C 

NP=NLIB+1 

INDICE = J(NLIB) 

M = J(NLIB)-1 

DO 1 I=2,NLIB 

IB=NP-I 

M = M + J(IB) 

IB=M+I 

INDICE = INDICE + IBIN(IB,I) 
1 CONTINUE 
RETURN 
END 

Subroutine ESPGN is the inverse of the indexing function. Given the index N it 
calculates the array J of dimension NLIB which contains the exponents. The value of 
N must be positive (not checked) and must not exceed the maximal index implicitly 
introduced by the initial choice of NLIB and NORD passed to BINOM. The latter error is 
actually checked (this does not increase the computation time). The code implements 
the recursive algorithm described in sect. 5.4, again using iteration. 

SUBROUTINE ESPON(N, J,NLIB) 

C 

C Compute the exponents J correponding to the 

C index N. 



C 



COMMON /BINTAB/ IBIN (NBNl , NBN2) , NPIUl , NMENl , NFAT , NBIN 
DIMENSION J (NLIB) 

NM=NLIB-1 

NP=NLIB+1 

DO 1 K=NLIB,NFAT 

IF (N.LT. IBIN (K, NLIB)) GO TO 2 

CONTINUE 

WRITE(6,1000) 

STOP 

NN = K-1 

M = N-IBIN(NN,NLIB) 
IF(NLIB-2) 8,6,7 
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7 DO 3 I = 2,NM 
L = NP-I 

DO 4 K=L,NFAT 
IF(M.LT.IBIN(K,L)) GO TO 5 

4 CONTINUE 

5 IB=NLIB-L 

J (IB) = NN-K 
NN = K-1 

M = M - IBIN(NN,L) 
3 CONTINUE 

6 J(NM) = NN-M-1 
J(NLIB) = M 
RETURN 

8 J(1)=N 
RETURN 

1000 FORMAT (//,5X,15HERR0R SUB ESPON,//) 

END 

The code described here is the skeleton of a program performing algebraic manipu- 
lation on polynomial. Such a program must include a call to BINOM in order to initialize 
the table of binomial coefficients. 

In order to store the coefficient of a monomial with exponents J (an integer array 
with dimension NLIB the user must include a statement like 

K = INDICE(J,NLIB) 

and then say, e.g., F(K)=. . . which stores the coefficient at the address K of the array F. 

Suppose instead that we must perform an operation on all coefficients of degree 
lORD of a given function F. We need to perform a loop on all the corresponding indices 
and retrieve the corresponding exponents. Here is a sketch of the code. 

C Compute the minimum and meiximum index NMIN and NMAX 

C of the coefficients of order lORD. 

C 

IB=I0RD+NMEN1 

NMIN = IBIN(IB,NLIB) 

IB=IORD+NLIB 

NMAX = IBIN(IB,NLIB) - 1 

C 

C Loop on all coefficients 

C 

DO 1 N = NMIN, NMAX 
CALL ESPON (N, J, NLIB) 

. . . more code to operate on the coefficient F (N) . . . 
1 CONTINUE 

Let us add a few words of explanation. According to (20), the index of the first coefficient 
of degree s in n variables is /(s, 0, . . . , 0) = N{n^ s — 1), and we also have A^(n, s — 1) = 
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("^n~^) explained at the end of sect. 5.3. This explains how the limits NMIN and 
NMAX are calculated as N{n, s — 1) and N{n, s + 1) — 1, respectively. The rest of the code 
is the loop that retrieves the exponents corresponding to the coefficient of index N. 



6. Trigonometric polynomials 

Let us now consider the more general case An = • The index n in An denotes again the 
dimension of the space. The name used in the title of the section is justified because this 
case occurs precisely in the representation of trigonometric polynomials, as explained 
in sect. 3.2. 

We shall now denote by = + . . . + the length (or norm) of the vector 
k e Z". The tail t{k) of a vector k will be defined again as t{k) = {k2, . . . , kn). 

6. 1 Ordering relation 

Pick a fixed n, and consider the finite family of sets = Z , . . . , An = Z"'. 

Let k, k' e Am, with any 1 < m < n. We shall say k' ~< k in case one of the following 
conditions is true: 



(i) 


m 


> 1 


A 


\k'\ 


< \k\ 








(ii) 


m 


> 1 


A 


\k'\ 


= \k\ 


A 


l^'il 


> \ki\ ; 


(in) 


m 


> 1 


A 


\k'\ 


= \k\ 


A 


1^11 


= \ki\ A k[> ki; 


(iv) 


m 


> 1 


A 


\k'\ 


= |A;| 


A 


K- 


= fci A t{k') -< t{k) 



In table 3 the order resulting from this definition is illustrated for the cases m = 2, 3, 4. 

If m = 1 this ordering coincides with the ordering in Z introduced in sect 4. For 
m > 1, if (i), (ii) and (iii) do not apply, then (iv) means that one must decrease the 
dimension n by replacing k with its tail t{k), and retry the comparison. Eventually, one 
ends up with m = 1, falling back to the one dimensional case to which only (i) and (iii) 
apply. 

The ordering in this section has been defined for the case An = Z"'. However, it 
will be useful to consider particular subsets of 7/^. The natural choice will be to use 
again the ordering relation defined here. For example, the case of integer vectors with 
non negative components discussed in sect. 5.1 can be considered as a particular case: 
the restriction of the ordering relation to that case gives exactly the order introduced in 
sect. 5.1. Just remark that the condition (iii) above becomes meaningless in that case, 
so that it can be removed. 

The set Vn{k) of the elements preceding k G A^ in the order above is defined as 
in sect. 5.1. Following the line of the discussion in that section it is also convenient to 
give some more definitions. Pick a vector k G An, and define the sets Bn\k), Bn^\k), 
Bn"\k) and Bn^\k) as the subsets of An satisfying (i), (ii), (iii) and (iv), respectively. 



Table 3. Ordering of integer vectors in Z'" for m — 2,3, 4. 
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( 0, 


3) 


( 0, 0, 2) 


( 0, 


2, 


0, 


0) 


24 


( 0: 


-3) 


( 0, 0,-2) 


( 0, 


-2, 


0, 


0) 



in the ordering algorithm above. FormaUy, 



Bli\k) = {k'eAn 
(25) B^^'Hk) = {k' e An 
Bi:^^\k) = {k' e An 
Bli'^Hk) ^ {k' e An 



\k'\ <\k\}, 

\k'\ = \k\ A \k[\ > \ki\} , 

\k'\ = \k\ A \k[\ = \ki\ A k[ > ki} , 

\k'\ ^ \k\ A k[^ki A t{k') <t{k)} . 
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The sets B^\k), Bli'\k), Bt'\k) and B^^''\k) are pair wise disjoint, and moreover 

B« {k) U Bf^ (k) U Bi"^) (k) U Bi-) (/c) = T'l/c) . 
This easily foUows from the definition. 

6.2 Indexing function for trigonometric polynomials 

Let k e An- In view of the definitions and of the properties above the index function, 
defined as in (14), turns out to be 

(26) /(O) = , I{k) = #B« (k) + i^B^^ (k) + #B4-^) (k) + #Bli^^ (k) . 

Let us introduce the J-table and the A'"-table as 
J{n,s) = #{A; e An : |A;| = s} , 

s 

N{n,s) = ^J{n,j) forn>l,s>0. 



(27) 



We claim that the index function can be recursively computed as 
(28) 
7(0) = , 

(N{n,\k\-1) for = A /ci > , 

N(n,\k\ + 1 for = j/cl A /ei < , 

N{n,\k\-l) + N{n-l,\k\-\ki\-l)+I{t{k)) for |fci| < |fc| A fci > , 

\.N{n,\k\-l) + N{n-l,\k\-\ki\) + I{t{k)) for |A;i| < |/c| A ki < . 



m 



This formula follows from 
(29) 

#Bf\k) = 



(30) 

(31) 
(32) 



#B«(fc) = iV(n,|/c|-l) ; 




#Bt'\k) = 
i^B^^\k)^ 







J{n-l,\k\-\k^\) 

jim) 



for 


l^il 




\k\ 


1 




1) for 


l^il 


< 


\k\ 






for 




< 


\k\ 


A ki 


>o. 


for 


\ki\ 


< 


\k\ 


A ki 


<0; 


for 


\ki\ 




\k\ 






-\ki\-l) for 




< 


\k\ 







The equality (29) is a straightforward consequence of the definition (25). The equal- 
ity (30) follows by remarking that for \ki \ — \k\ we have Bn^\k) — 0, and for < 
we have 

= B+{k) U B-{k) , S+(fc) n i?-(/c) = , 
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with 

B+{k) = U {k'^An : K ^\k\-l A \t{k)\ = /} , 

0<Z<|fc|-|fci| 

B-{k)^ \J {k' e An : k[ ^ I - \k\ A \t{k)\ ^ 1} ; 

0<l<\k\-\ki\ 

use also #B^{k) — #B~{k). The equahty (31) foUows from 
Br)(/c) = 



l^ll 


= \k\ 


A ki 


>o , 




= \k\ 


A ki 


< . 



Coming to (32), remark that 

Bli'^Hk) = {k' e An : \tik')\ = \k\ -\ki\ A t{k') ■< t{k)} . 
Proceeding as in the polynomial case we find again 

#,Bi-)(/c) = #{A e An-i : |A| = \k\-\ki\ A t{k)} , 
and (32) follows by remarking that 

Vn-i{t{k)) = {A e An-i ■■ |A| = |A;| - \ki\ A A -< t{k)}U{\ E An-i ■ |A| < \k\ - \ki\} . 
Adding up all contributions (28) follows. 

6.3 Construction of the tables 

We show now how to construct recursively the J-table, so that the A'"-table can be 
constructed, too. For n = 1 we have, trivially, J(l, 0) = 1 and J(l, s) = 2 for s > 0. For 
n > 1 use the elementary property 

{keAn ■ \k\ =s}= y {keAn-ki=jA \t{k)\ = s- \j\} . 

—s<j<s 



Therefore 



(33) 



J(1,0) = 1 , 
J(l,s) = 2 , 

s 

J{n, s) — ^ J(n — l,s — for n > 1 . 



This completely determines the J-table. 
6.4 Inversion of the index function 

The problem is to find the vector k of given dimension n corresponding to the given 
index I. For n = 1 the function I{k) and its inverse are given by (15) and (16). 

Therefore in the rest of this section we shall assume n > 1. We shall give a recursive 
algorithm, showing how to determine ki and l(t{k)). 
(i) If ^ = then k — 0, and there is nothing else to do. 
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(ii) Assuming that Z > 0, determine s such that 

N{n,s- 1) <l < N{n,s) . 

From this we know that = s. 

(iii) Define h=l- N{n, s - 1), so that l{t{k)) < h by (28). Ii h = set si = 0; else, 
determine s' such that 

N(n -l,s' -1) <h < N(n - 1, s') , 

and let si = min(s',s). In view of l(t{k)) < h we know that \t{k)\ < si. Remark 
also that si = if and only if li = 0. For, if si > 1 then we have h > N{n—1, 0) = 1. 

(iv) If h = 0, then by the first of (28) we conclude 

ki = \k\=s , t{k) = , 

and there is nothing else to do. 

(v) If = 1, then by the second of (28) we conclude 

ki = —\k\ = —s , t{k) = , 

and there is nothing else to do. 

(vi) If Zi > 1 and si > 0, we first look if we can set < fci < In view of the third 
of (28) we should have 

\k\-ki=si, \t{k)\^si, I{t{k)) =li- N{n-l,si-l) . 

This can be consistently made provided the conditions 

si > and l{t{k)) > N{n - 1, si - 1) 

are fulfilled. The condition s > is already satisfied. By (28), the second condition 
is fulfilled provided li > 2N{n — 1, Si — 1). This has to be checked, 
(vi.a) If the second condition is true, then set /ci = | A; | — s'l, and recall that |t(/c)| — Si. 
Hence, we can replace n, I, and s by n — 1, li — N{n — 1, si — 1) and si, 
respectively, and proceed by recursion restarting again from the point (iii). 
(vi.b) If the second condition is false, then we proceed with the next point. 

(vii) Recall that h > 1, and remark that we have also Si > 1. Indeed, we already know 
si > 0, so we have to exclude the case si — 1. Let, by contradiction, si = 1. Then 
we have h >2 = 2N{n — 1, si — 1), which is the case already excluded by (vi). We 
conclude si > 1. We look now for the possibility of setting < and ki < 0. 
In view of the fourth of (28) we should have 

\k\ + ki = si-l, \t{k)\=si-l, l{t{k)) =h- N{n-l,si-l) . 

This can be consistently made provided the conditions 

si > 1 and I{t{k)) > N{n-l,si-2) 

are fulfilled. The condition si > 1 is already satisfied. As to the second condition, 
by (28) it is fulfilled provided h > N{n - 1, si - 1) + N{n - 1, si - 2). This has to 
be checked. 
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(vii.a) If the second condition is true, then set /ci = — + si — 1, and recaU that 
\t{k) I — si — 1. Hence, we can replace n, /, and s by n — 1, /i — N{n — 1, si — 2) 
and si — 1, respectively, and proceed by recursion restarting again from the 
point (iii). 

(vii.b) If the second condition is false we must decrease Si by one and start again with 
the point (vi); remark that Si > 1 implies Si — 1 > 0, which is the first of the 
two conditions to be satisfied at the point (vi), hence the recursion is correct. 
Since Zi > 1 we have h > 2N{n—l, 0), so that the conditions of point (vi) are satisfied for 
s = 1. Hence the algorithm above does not fall into an infinite loop between points (vi) 
and (vii). On the other hand, for n = 1 either (iii) or (iv) applies, so that the algorithm 
stops at that point. 



7. Storing the coefficients for sparse functions 

The method of storing the coefficient using the index, as illustrated in sect. 3, is the 
most direct one, but reveals to be ineffective when most of the coefficients of a function 
are zero (sparse function). For, allocating memory space for all coefficients results in a 
wasting of memory. 

A method that we often use is to store the coefficients using a tree structure based 
on the index. However we should warn the reader that the method described here has 
the advantage of being easily programmed, but does not pretend to be the most effective 
one. Efficient programming of tree structure is described, e.g., in the monumental books 
The art of computing programming, by D.E. Knuth [20]. 

7.1 The tree structure 

The first information we need is how many bits are needed in order to represent the 
maximum index for a function. We shall refer to this number as the length of the index. 
In the scheme that we are presenting here this is actually the length of the path from 
the root of the tree to its leave, where the coefficient is found. 

In fig. 1 we illustrate the scheme assuming that 4 bits are enough, i.e., there are 
at most 16 coefficients indexed from to 15. The case is elementary, of course, but the 
method is the general one, and is extended to, e.g., several millions of coefficients (with 
a length a little more than 20) in a straightforward manner. The bits are labeled by 
their position, starting from the less significant one (choosing the most significant one 
as the first bit is not forbidden, of course, and sometimes may be convenient). The label 
of the bit corresponds to a level in the tree structure, level being the root and level 
3 being the last one, in our case. At level zero we find a cell containing two pointers, 
corresponding to the digit and 1, respectively. To each digit we associate a cell of 
level 1, which contains a pair of pointers, and so on until we reach the last level (3 in 
our case). Every number that may be represented with 4 bits generates a unique path 
along the tree, and the last cell contains pointers to the coefficient. The example in the 
figure represents the path associated with the binary index 1010, namely 10 in decimal 
notation. 
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Figure 1. Illustrating the tree structure for a 4-bit long index (see text). 



Let us also illustrate how this structure may be represented in memory, trying to 
avoid wasting of space. We use two separate arrays, the first one for pointers and the 
second one for the coefficients, as illustrated in fig. 2. The cells containing pairs of 
pointers are allocated in the first array, the root of the tree having label zero. The label 
of a cell is always even: the first element corresponds to the zero bit, the next one (with 
odd label) to the bit one. 

The arrays are initially allocated with appropriate size, and are cleared. A good 
method is to fill the array of pointers with —1 (denoting an unused pointer) and the 
coefficients table with zeros. We also keep track of the first unused cell in the array, 
which initially is set to 2 because the root cell is considered to be in use, and of the first 
free coefficient, which initially is 0. 

We shall use the following notations: cell(2j) is the cell with even label 2j in the 
array; cell(2j, 0) and cell(2_;/, 1) are the pointers corresponding to a bit or 1 which 
are stored at locations 2j and 2j + 1, respectively, in the array of pointers; coef (j) is 
the j-th element of the array of coefficients; cc is the current cell and cb is the current 
bit (see below for the meaning); fp is the label of the first free (unused) cell of pointers; 
f c is the label of the first free coefficient; £ is the length of the index. 

7.2 Storing the Brst coeiEcient 

Let us describe how the first coefficient is stored. Suppose we want to store the value x 
as the coefficient corresponding to a given index. Here is the scheme. 

(i) Initialization: set cc = and cb = 0. The values of f p = 2 and f c = have already 
been set when during the array allocation. 

(ii) Creating a path: repeat the following steps until cb equals £ — 1: 
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Figure 2. Illustrating how the tree structure is stored in memory (see text). 



(ii.a) if the bit at position cb in the index is 0, then redefine cell(cc, 0) = f p; else 
redefine cell(cc, 1) = fp; 

(ii.b) set cc = fp and increment fp by 2 (point to the next free cell); 

(ii.c) increment cb by 1 (next bit). 

(iii) Store the coefRcient: 

(iii.a) if the bit at position cb in the index is 0, then redefine cell(cc,0) = fc; else 
redefine cell(cc, 1) = f c; 

(iii.b) set coef (f c) = x; 

(iii.c) increment f c by 1 (point to the next free coefficient). 

Programming this algorithm in a language such as C or FORTRAN requires some 10 to 20 
statements. 



Let us see in detail what happens if we want to store the coefficient 0.6180339 with 
index 1010 and £ — 4, as illustrated in fig. 2. Here is the sequence of operations actually 
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made 



step (i): 


cc = , 




cb = , 


fp = 2 , 


f c = 


0; 




step (ii): 


cell(0, 0) 


= 2 , 


cc = 2 , 


fp = 4 , 


cb = 


1 , 


then , 




cell(2, 1) 


= 4, 


cc = 4 , 


fp = 6 , 


cb = 


2 , 


then , 




cell(4, 0) 


= 6 , 


cc = 6 , 


fp = 8 , 


cb 


3, 


end of loop ; 


step (iii): 


cell(6, 1) 


= , 


coef(0)= 


=0.6180339 , 


fc = 


1 , 


end of game 



After this, the contents of the arrays are as represented in fig. 2. 
7.3 Retrieving a coeiRcient 

The second main operation is to retrieve a coefficient, which possibly has never been 
stored. In the latter case, we assume that the wanted coefficient is zero. Here is a scheme. 

(i) Initialization: set cc = and cb = 0. 

(ii) Follow a path: repeat the following steps until cb equals £: 
(ii.a) save the current value of cc; 

(ii.b) if the bit at position cb in the index is 0, then redefine cc as cell(cc,0); else 

redefine cc as cell(cc, 1); 
(ii.c) if cc = — 1 then the coefficient is undefined. Return as the value of the 

coefficient; 
(ii.d) increment cb by 1 (next bit). 

(iii) CoefEcient found: return the coefficient coef (cc). 

Let us give a couple of examples in order to better illustrate the algorithm. Suppose 
that we are looking for the coefficient corresponding to the binary index 1010. By 
following the algorithm step by step, and recalling that in our example the length of 
the index is 4, the reader should be able to check that the sequence of operations is the 
following: 

step (i): cc = , cb = ; 

step (ii): cc = 2 , cb = 1 , then , 

cc = 4 , cb = 2 , then , 

cc = 6 , cb = 3 , then , 

cc = cb = 4 , end of path ; 
step (iii): return 0.6180339 . 

The returned value is that of coef (0), stored in the location of the coefficients array. 

Suppose now that we are looking for the coefficient corresponding to the binary 
index 1110. Here is the actual sequence of operations: 



step (i): 


cc 


= 0, 


cb = 


0; 




step (ii): 


cc 


= 2, 


cb = 


1 , 


then , 




cc 


= 4, 


cb = 


2, 


then , 




cc 


= -1 , 


cb = 


2, 


return zero 



Here the algorithm stops because a coefficient has not been found. 
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Figure 3. Inserting a new coefficient in a tree structure (see text). 

7.4 Other operations 

Having implemented the two operations above, the reader should be able to implement 
also the following operations: 

(i) storing a new coefficient corresponding to a given index; 

(ii) adding something to a given coefficient; 

(iii) multiplying a given coefficient by a number. 

These are the basic operations that we need in order to perform an elementary computer 
algebra. Let us add a few hints. 

Storing a new coefficient requires perhaps some moment of thinking. Using the 
index, one should follow the corresponding path in the tree (as in the operation of 
retrieving) until either happens: the coefficient is found, or the search fails at some 
point. If the coefficient is found, then it can be overwritten if the new value has to 
replace the old one. On failure, the path must be completed by appropriately defining 
the pointers (as in the case of the first coefficient), and then the coefficient can be stored 
in the appropriate location. As an exercise, suppose that we want to store the coefficient 
1.4142136 corresponding to the binary index 1110. After completing the operation the 
memory should look as in fig. 3. 

Adding something to a given coefficient is not very different from the previous 
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operation. Just follow the path. If the coefficient is found, then add the wanted value 
to it. On failure, just change the "add" operation to a "store" one, and proceed as in 
the case (i). 

Multiplying a coefficient by a constant is even easier. If the coefficient is found, then 
do the multiplication. On failure, just do nothing. 

Further operations can be imagined, but we think that we have described the basic 
ones. There are just a couple of remarks. 

The method illustrated here uses an amount of memory that clearly depends on 
the number of non zero coefficients of a function. However, this amount is typically not 
known in advance. Thus, enough memory should be allocated at the beginning in order 
to assure that there is enough room. When a function is filled, and we know that it 
will not be changed, the excess of memory can be freed and reused for other purposes. 
Every operating system and language provides functions that allow the programmer to 
allocate memory blocks and resize them on need. 

A second remark is that other storing methods can be imagined. E.g., once a func- 
tion is entirely defined it may be more convenient to represent it as a sequential list of 
pairs (index, coefficient). This is definitely a very compact representation for a sparse 
function (although not the best for a crowded one). 

8. Applications 

We report here some examples of application of algebraic manipulation that have been 
obtained by implementing the formal algorithm of sect. 2. We consider three cases, 
namely the model of Henon and Heiles, the Lagrangian triangular equilibria for the Sun- 
Jupiter system and the planetary problem including Sun, Jupiter, Saturn and Uranus 
(SJSU). 

8. 1 The model of Henon and Heiles 

A wide class of canonical system with Hamiltonian of the form 

(34) H{x, y) = Y(y? + xl) + ^(yi + xl) + xjx^ 

has been studied by Contopoulos, starting at the end of the fifties, for different values 
of the frequencies. This approximates the motion of a star in a galaxy, at different 
distances from the center. A wide discussion on the use of these models in galactic 
dynamics and on the construction of the so called "third integral" can be found in the 
book [6]. The third integral is constructed as a power series $ = $2 + ^3 + - • • where is 
a homogeneous polynomial of degree s which is the solution of the equation {H, = 0, 
where {•, •} is the Poisson bracket (see, e.g., [27] or [5]). A different method is based on 
the construction of the Birkhoff normal form [1]. 

A particular case with two equal frequencies and Hamiltonian 

(35) H{x, y) ^\{yl + xl) + \{yl + xl) + xix2 - ^xl 
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Figure 4. Poincare sections for the Henon and Heiles model. The energies are as 
in the original paper. 



has been studied by Henon and Heiles in 1964 [15]. This work has become famous since 
for the first time the existence of a chaotic behavior in a very simple system has been 
stressed, showing some figures. It should be remarked that the existence of chaos had 
been discovered by Poincare in his memory on the problem of three bodies [22], but it 
had been essentially forgotten. 

A program for the construction of the third integral has been implemented by 
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Figure 5. Level lines of the first integral truncated at orders 4 and 8, for energy 
E = j^. The figure for truncation orders up to 58 are actually the same as for 
order 8. 

Contopoulos since 1960. He made several comparisons between the level lines of the 
integral so found on the surface of constant energy and the figures given by the Poincare 
sections of the orbits. A similar calculation for the case of Henon and Heiles has been 
made by Gustavson [14], who used the normal form method. The third integral was 
expanded up to order 8, which may seem quite low today, but it was really diflBcult to 
do better with the computers available at that time. Here we reproduce the figures of 
Gustavson extending the calculation up to order 58, which is now easily reached even 
on a PC. 

In fig. 4 we show the Poincare sections for the values of energy used by Henon and 
Heiles in their paper. As stressed by the authors, an essentially ordered motion is found 
for £■ < y^, while the chaotic orbits become predominant at higher energies. 

The comparison with the level lines of the third integral at energy i? = ^ is 
reported in fig. 5. The correspondence with the Poincare sections is evident even at 
order 8, as calculated also by Gustavson. We do not produce the figures for higher 
orders because they are actually identical with the one for order 8. This may raise the 
hope that the series for the first integral is a convergent one. 

Actually, a theorem of Siegel states that for the Birkhoff normal form divergence is 
a typical case [26] . A detailed numerical study has been made in [7] and [8] , showing the 
mechanism of divergence. Moreover, it was understood by Poincare that perturbations 
series typically have an asymptotic character (sec [23], Vol. II). Estimates of this type 
have been given, e.g., in [11] and [12]. 

For energy E = ^ (fig. 6) the asymptotic character of the series starts to appear. 
Indeed already at order 8 we have a good correspondence between the level lines and the 
Poincare section, as was shown also Gustavson's paper. If we increase the approximation 
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Figure 6. Level lines of the first integral truncated at orders 8, 32, 43 and 58, 

for energy E = A good correspondence with the Poincare sections is found at 
orders, roughly, 8 to 32. Then the level lines start to disprove, in agreement with 
the asymptotic character of the series. 

we see that the correspondence remains good up to order 32, but then the divergence 
of the series shows up, since at order 43 an unwanted "island" appears on the right side 
of the figure which has no correspondent in the actual orbits, and at order 58 a bizarre 
behavior shows up. 

The phenomenon is much more evident for energy E = ^ (fig. 7). Here some rough 
correspondence is found around order 9, but then the bizarre behavior of the previous 
case definitely appears already at order 27. 
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Figure 7. Level lines of the first integral truncated at orders 8, 9, 10 and 27, for 
energy E = ^. Some correspondence with the Poincare sections is found around 
the order 9. Then the level lines are definitely worse, making even more evident 
the asymptotic character of the series. 

The non convergence of the normal form is illustrated in fig. 8. Writing the homo- 
geneous terms of degree s of the third integral as $s = ^ {pj^kX-'y'^ , we may introduce 
the norm 

11*^11 = Zli'^j'fci • 

Then an indication of the convergence radius may be found by calculating one of the 
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Figure 8. The convergence radius evaluated with the ratio (left) and the root 
(right) criterion. In both cases the non convergence of the series is evident. 



quantities 



,1/s 




The first quantity corresponds to the root criterion for power series. The second one 
corresponds to the ratio criterion. The third one is similar to the ratio criterion, but in 
the present case turns out to be more effective because it takes into account the peculiar 
behavior of the series for odd and even degrees. The values given by the root criterion 
are plotted in the left panel of fig. 8. The data for the ratio criterion are plotted in the 
right panel, where open dots and solid dots refer to the second and third quantities in 
the formula above, respectively. In all cases it is evident that the values steadily increase, 
with no tendency to a definite limit. The almost linear increase is consistent with the 
behavior ||^s|| ~ s! predicted by the theory. 

8.2 The Trojan asteroids 

The asymptotic behavior of the series lies at the basis of Nekhoroshev theory on expo- 
nential stability. The general result, referring for simplicity to the case above, is that in 
a ball of radius g and center at the origin one has 

|$(t) - $(0)| < 0{g^) for \t\ < 0(exp(l/^")) , 

for some positive a < 1. This is indeed the result given by the theory (see, e.g., [11]). 
In rough terms the idea is the following. Due to the estimate ~ s\ and remarking 
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Figure 9. The estimated stability time and the optimal truncation order for the 
L4 point of the Sun-Jupiter system. 



that $ = {i/, $} starts with terms of degree s + 1, one gets |$| = 0{s\g'+^). Then 
one looks for an optimal degree s which minimizes the time derivative, i.e., s ^ 1/g. By 
truncating the integrals at the optimal order one finds the exponential estimate. 

However, the theoretical estimates usually give a value of g which is useless in prac- 
tical applications, being definitely too small. Realistic results may be obtained instead if 
the construction of first integrals for a given system if performed by computer algebra. 
That is, one constructs the expansion of the first integral up to an high order, compat- 
ibly with the computer resources available, and then looks for the optimal truncation 
order by numerical evaluation of the norms. 

The numerical optimization has been performed for the expansion of the Hamil- 
tonian in a neighborhood of the Lagrangian point L4, in the framework of the planar 
circular restricted problem of three bodies in the Sun-Jupiter case. This has a direct 
application to the dynamics of the Trojan asteroids (see [13]). 

The two first integrals which are perturbations of the harmonic actions have been 
constructed up to order 34 (close to the best possible with the computers at that time). 
The estimate of the time of stability is reported in fig. 9. The lower panel gives the 
optimal truncation order vs. log^g g- In the upper panel we calculate the stability time 
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Figure 10. The estimated stability time and the optimal truncation order for 
the SJSU planar system. The dashed line corresponds to the estimated age of the 
Universe. 

as follows: for an initial datum inside a ball of radius qq we determine the minimal 
time required for the distance to increase up to 2go. Remark that the vertical scale is 
logarithmic. The units are chosen so that ^ = 1 is the distance of Jupiter from the 
Sun, and t = 27r is the period of Jupiter. With this time unit the estimated age of the 
universe is about 10^. The figure shows that the obtained data are already realistic, 
although, due to the unavoidable approximations, only four of the asteroids close to L4 
known at the time of that work did fall inside the region of stability for a time as long 
as the age of the Universe. 
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8.3 The SJSU system 

As a third application we consider the problem of stability for the planar secular plan- 
etary model including the Sun and three planets, namely Jupiter, Saturn and Uranus. 
The aim is evaluate how long the semi-major axes and the eccentricities of the orbits 
remain close to the current value (see [25]). 

The problem here is much more difficult than in the previous cases. The Hamiltonian 
must be expanded in Poincare variables, and is expressed in action-angle variables for the 
fast motions and in Cartesian variables for the slow motions, for a total of 9 polynomial 
and 3 trigonometric variables. The expansion of the Hamiltonian in these variables 
clearly is a major task, that has been handled via computer algebra. 

The reduction to the secular problem actually removes the fast motions, so that 
we get an equilibrium corresponding to an orbit of eccentricity zero close to a circular 
Keplerian one, and a Hamiltonian expanded in the neighborhood of the equilibrium, 
which is still represented as a system of perturbed harmonic oscillators, as in the cases 
above. Thus, after a long preparatory work, we find a problem similar to the previous 
one, that can be handled with the same methods. 

The results are represented in fig. 10, where we report again the optimal truncation 
order and the estimated stability time, in the same sense as above. The time unit here is 
the year, and the distance is chosen so that = 1 corresponds to the actual eccentricity 
of the three planets. The result is still realistic, although a stability for a time of the 
order of the age of Universe holds only inside a radius corresponding roughly to 70% of 
the real one. 
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